Absence of Thermodynamic Phase Transition in a Model Glass Former 
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The glass transition can simply be viewed as the point at which the viscosity of a structurally 
disordered liquid reaches 10 13 Poise Q. This definition is operational but it sidesteps fundamental 
controversies about the glass: Is the transition a purely dynamical phenomenon pi? This would 
mean that ergodicity gets broken, but that thermodynamic properties of the liquid remain un- 
changed across the transition if determined as thermodynamic equilibrium averages over the whole 
phase space. The opposite view claims that an underlying thermodynamic phase transition is 
responsible for the dramatic slowdown at the liquid-glass boundary. Such a phase transition (which 
shows up in proper equilibrium phase space averages) would trigger the dynamic standstill, and 
then get masked by it. A recent Monte Carlo algorithm Q introduces non-local moves for hard-core 
systems in a way which preserves micro-reversibility. Here we use this method to equilibrate a 
two-dimensional hard disc system far within its glassy phase. We show that indications of a ther- 
modynamic transition are lacking up to very high densities, and that the glass is indistinguishable 
from the liquid on purely thermodynamic grounds. 



We considered polydisperse hard discs i = 1,...,N 
with radii (r^ — fj_i = A/(iV — 1)) in a square box of 
volume V — L 2 with periodic boundary conditions. The 
polydispersity A/ri was kept fixed and the system stud- 
ied as a function of the density p — nJ2i r i/L 2 - Note 
that hard-core systems are athermal: the phase diagram 
and the dynamics (up to a trivial rescaling) are indepen- 
dent of temperature. The external control parameters 
are the density (for the [NV] ensemble) or the pressure 
(for [NP]). 

Conventional, local-move Monte Carlo simulations in 
the [NV] ensemble were used to compute effective diffu- 
sion constants || at densities for which the local algo- 
rithm was still ergodic, cf. fig. 1. We failed to detect any 
finite-size effects by comparing runs with 256 and 1024 
discs at density p = 0.764, but we noticed strong depen- 
dence of the diffusivities on the disc size. For each disc i 
we found its diffusivity Di to agree very well with [fl 
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For the ~ 180 largest discs we found no systematic de- 
pendence of pf = p G ^ 0.805 ± 0.01 (with a best fit 
a ~ 2.4). For the smallest discs i, the extrapolated val- 
ues for pf increased slowly. Our result for p agrees with 
what was found in a related system (B. Doliwa, Diploma 
Thesis, University of Mainz (1999)). 

We next performed simulations with the cluster Monte 
Carlo algorithm [Q. There, groups of discs are swapped 
around a 'pivot' (cf. fig. 4) in a way which represents 
an alternative Markov-chain sampling of the Boltzmann 
distribution. The correctness of both implementations 
was validated by comparing structural quantities (pair 
correlation functions) and thermodynamic variables in 
the manifestly liquid regime, where the local algorithm 
still converges well. We found that the cluster Monte 
Carlo algorithm does not slow down as we pass pa- 



FIG. 1. Effective diffusion constants as a function of the 
radius for polydisperse hard discs. Upper 6 curves: 256 
discs with polydispersity A/ri = 19 were simulated using 
the local-move Monte Carlo algorithm at densities at which 
the algorithm equilibrated. The effective diffusion constants 
were calculated via Di =< (xi(to + t) — Xi(to)) 2 > /(4t), for 
|xi(to + t) ~ x i(^o)| S> Hv- The time was measured in Monte 
Carlo sweeps. Only the motion relative to the largest disc was 
taken into account. These asymptotic constants agree (up to 
a rescaling of the time unit) with what would be obtained 
in a molecular dynamics simulation [6]. Lowest curve: as 
above, but equilibrated initial configurations were provided 
by the cluster Monte Carlo algorithm. We observe that the 
smallest ~ 50 discs are able to pass through the blocked ma- 
trix made up of the largest ~ 200 discs. These findings are 
well consistent with the estimated values of the transition 
densities. 
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FIG. 2. Equilibrium compressibility. System as before was 
simulated using the cluster Monte Carlo algorithm. Equilib- 
rium equation of state (inset) and inverse compressibility ft -1 
(main figure) were obtained in the [NP] ensemble, k, was ob- 
tained from the fluctuations of the volume (points with error 
bars) and from the derivative of the best V(P) fit according 
to the functional form V(P) — const — aP~ b (Parameters of 
the best fit are const — 84.3, a = 254, b = 0.69). Indications 
of an equilibrium phase transition are lacking. The gray lines 
indicate the localization of the glass transition as extrapolated 
from fig. 1. 

We used the cluster Monte Carlo algorithm to com- 
pute the equation of state in the [NP] ensemble S . For 
pressures corresponding to the glass transition density 
p G , we reached a relative precision for V of .05% during 
a one-week simulation on a single processor workstation 
(JV = 256). The compressibility of the system was ob- 
tained both by deriving the algebraic fit to V(P) (cf. 
inset of fig. 2) and by computing the volume fluctua- 
tions: 

-1 d(V) _ (v 2 ) - (V) 2 
- (V) dP ~ (V) [ ' 

The agreement of the two methods (main fig. 2) tes- 
tifies to the excellent convergence of the cluster Monte 
Carlo algorithm. Some of the compressibility data 
at high density were obtained from several indepen- 
dent runs with identical results. In contrast, the local 
Monte Carlo algorithm yielded non-reproducible (gener- 
ally lower) compressibilities. 

The results presented in fig. 2 allow us to draw far- 
reaching conclusions: 

i) In the polydisperse hard-disc system considered, the 
glass is indistinguishable from the liquid on purely ther- 
modynamic grounds, as we find no indication of a phase 
transition. 

ii) Furthermore, the glass remains thermodynamically 
stable. No tendency towards crystallization or phase sep- 
aration was detected. At smaller values of the polydis- 
persity A/n, a crystalline state is present, and can be 



readily detected in the equation of state jlOj. In a re- 
lated two-dimensional system of binary mixtures, phase 
separation was evidenced by the cluster Monte Carlo al- 
gorithm 

Hi) Finally, the smooth behaviour of the compressibil- 
ity as we cross p G indicates that the ergodicity of the clus- 
ter Monte Carlo algorithm is preserved up to the highest 
densities studied. 
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FIG. 3. Distribution of inherent structure densities for 15 
discs with polydispersity A/ri = 19. Equilibrium configura- 
tions were generated during a long run of the cluster Monte 
Carlo algorithm. The corresponding inherent structures were 
obtained by iterating a rescaling of the particle radii and lo- 
cal Monte Carlo moves ( "rattling" ) until a jamming state (of 
density pis) was reached. As p increases, inherent structures 
with higher pis become more probable. The same inherent 
structures with very large values of pis were found in the 
simulation at p = 0.82 and p = 0.84. 

Beyond the above point Hi), the cluster Monte Carlo 
algorithm has passed an extremely stringent ergodicity 
test in a system with 15 polydisperse discs (where we 
also identified a density p of complete dynamical stand- 
still) : A large number of equilibrated configurations were 
stored during single cluster-simulation runs in the [NV] 
ensemble at various densities p. For each stored config- 
uration we then alternated very small up-scaling of each 
radius r,; — > (1 + e)ri (i = 1,...,N) with a few lo- 
cal Monte Carlo moves ("rattling"). The procedure was 
repeated until a jamming state (with unchanged poly- 
dispersity A/ri) was reached. This state represents an 
inherent structure (IS), generalised from what is done in 
thermal systems by a quench to zero temperature |l2],[l3| . 
Histograms of the IS densities pis for different values of 
p show that the algorithm explores the remaining regions 
of phase space as the density p is increased, cf. fig. 3. 
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FIG. 4. Example of a cluster move at high density. The 
(+) denote the randomly chosen 'pivot' with respect to which 
the cluster (grey discs) can be nipped [7]. Periodic boundary 
conditions are indicated. The inherent structure related to 
the configuration a) is conjectured to yield the closest packed 
state for N = 15, A/n = 19. 

A non-ergodic algorithm would get stuck in one or a 
few IS (during a single run) as regions of phase space 
become dynamically inaccessible. We observed the con- 
trary: almost identical sets of IS with extremely large 
values of pjs 3> p appeared in simulations at different, 
very high, values of p (for intermediate values of pis, the 
match could no longer be achieved, simply because of 
the very large number of different IS, even in a system 
with 15 discs, cf. fig. 3). Figure 4a) shows a config- 
uration at p = 0.86, which can be blown up into the 
densest IS which we were able to find at polydispersity 
A/n = 19 {pis = 0.8938). This IS (as many others) 
recurred frequently during a single run both at p = 0.84 
and at p = 0.82. We conjecture that the IS related to fig. 
4a) solves the closest-packing problem for N — 15 and 
A/n = 19. 

In conclusion, our results show that the new Monte 
Carlo algorithm (which can be generalised to incorporate 
a smooth potential [Q) works extremely well at densities 
at which the local method is totally stuck. Besides the 
calculation of thermodynamic quantities, the algorithm 
can also be used to study the finite-time dynamics (start- 
ing from equilibrated samples) at p ~ p G . This is how the 
lowest curve in fig. 1 was obtained. As the initial condi- 
tion is thermalised, we gain access to unbiased dynamical 
information in a regime which used to be out of reach of 
rigorous simulations because not all discs can be prop- 
erly equilibrated with the local algorithm. Clearly, the 
method represents a very powerful tool to study glasses. 
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